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Abstract 

We examine two models for hepatitis C viral (HCV) dynamics, one for monotherapy with interferon 
(IFN) and the other for combination therapy with IFN and ribavirin. Optimal therapy for both the models 
is determined using the steepest gradient method, by defining an objective functional which minimizes 
the infected hepatocyte levels, virion population and the side-effects of the drug(s). The optimal therapy 
for both the models shows an initial period of high efficacy, followed by a gradual decline. The period of 
high efficacy coincides with a significant decrease in the infected hepatocyte levels as well as viral load, 
whereas the efficacy drops after liver regeneration through restored hepatocyte levels. The period of high 
efficacy is not altered significantly when the cost coefficients are varied, as long as the side effects are 
relatively small. This suggests a higher dependence of the optimal therapy on the model parameters in 
case of drugs with minimal side effects. 

We use the Latin hypercube sampling technique to randomly generate a large number of patient 
scenarios (i.e, model parameter sets) and study the dynamics of each set under the optimal therapy 
already determined. Results show an increase in the percentage of responders (as indicated by drop in 
viral load below detection levels) in case of combination therapy as compared to monotherapy. Statistical 
tests performed to study the correlations between sample parameters and the time required for the viral 
load to fall below detection level, show a strong monotonic correlation with the death rate of infected 
hepatocytes, identifying it to be an important factor in deciding individual drug regimens. 

Keywords: Hepatitis C, Optimal Control, Steepest-Gradient Method, Latin Hypercube Sampling, 
Statistical Tests 

1 Introduction 

Hepatitis C (HCV) is an infectious disease which spreads through blood contact. It is estimated that 
HCV has infected 170 million individuals worldwide HI. Roe and Hall |2| estimate that about 50 — 80% 
of HCV infected cases are chronic in nature. Of these chronic cases, about 10 — 20% develop into liver 
cirrhosis of which, about 5% develop hepatocellular carcinoma (HCC). The extend of prevalence of HCV 
varies widely across geographical locations While the dominant mode of HCV transmission in United 

States, Europe and Australia is injecting drug use lO], the absence of reliable screening for HCV amongst 
blood donors remains a major challenge in combating the spread of the disease in India ||5]|6l. 

Mathematical modelling and quantitative analysis of hepatitis C infections has been explored extensively 
over the last decade. Most of the modelling has been restricted to the short term dynamics of the model. One 
of the earliest models was proposed by Neumann et al. Q, who examine the dynamics of HCV in presence 
of Interferon-a (IFN-a) treatment. They find that the primary role of IFN is in blocking the production 
of virions from the infected hepatocytes. However, IFN has little impact when it comes to controlling the 
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infection of the hepatocytes. Dixit et al. HI improved upon Q by including the effects of ribavirin , which 
in turn results in a fraction of the virions being rendered noninfectious. Their model is able to explain 
clinically observed biphasic decline patterns amongst patient population. Their study also shows that while 
IFN plays a pivotal role in the first phase decline of viral load, ribavirin has very little impact. However, 
in case of low IFN efficacy, ribavirin makes a significant contribution to the second phase of decline. The 
model could not successfully explain the triphasic decline patterns, as well as some cases of non-responders. 
Dahari et al. (91 in a subsequent and improved model, take into account the homeostatic mechanisms for the 
liver by incorporating a growth function. This model successfully explains the triphasic decline, as well as 
therapeutic failures. 

Control theory has found wide ranging applications in biological and ecological problems ifTOl . In 
biomedical problems, techniques from control theory are of great use in developing optimal therapeutic 
strategies. The treatment regimen is usually taken to be the control variable, with the aim of minimizing the 
detrimental effects of the medical condition. For instance, in cancer modeling, DePillis et al. ifTTI examine 
a mathematical model of tumor-immune interactions and determine the optimal chemotherapy with the goal 
of minimizing the tumor density, as well as the therapeutic side effects, for both linear and quadratic con- 
trol. Fister and Panetta |[T2l consider three different cell-kill models and determine the optimal drug dosage 
which minimizes the cancer mass and the cost of treatment for all three models. Murray [13] uses linear 
control approach for cancer chemotherapy models with toxicity limit. Optimal control theory can be applied 
to epidemic models, such as the SEIR model, with the goal of minimizing the number of infectious individ- 
uals and the overall cost of the vaccination programme lITOl . In case of treatment of diabetes, the control is 
taken as the insulin injection levels with the objective of minimizing the difference between current and the 
desired glucose levels. Chavez et al. |[T4l present an optimal insulin delivery for type 1 diabetic patients. 

Control theoretic approach has been applied extensively in case of virological models, especially in case 
of HIV models. One of the earliest papers in this area by Kirschner et al. |[T5l uses an existing ODE 
model which incorporates the dynamics of the immune system and HIV. Using an objective functional 
which maximizes the levels of T-cells and minimizes the cost of treatment, the optimal chemotherapy (like 
protease inhibitors) is determined for various stages of initiation of treatment. Joshi |[T6l considers a model 
with combination therapy and numerically solves the optimal system for the optimal treatment regimen. 
Adams et al. IITtI present the optimal treatment protocols (along with modeling and data analysis) for HIV 
dynamics. Stengel lITSl . in his paper, obtains the optimal therapies for drug resistant strains of HIV which 
evolve rapidly, as a result of high viral turnover and mutation rates. He observes that in such cases continued 
treatment is imperative for sustained remission. 

In the case of HCV, Chakrabarty and Joshi im consider a model (motivated by Q El Hi) for HCV 
dynamics under combination therapy of interferon and ribavirin. An objective functional is formulated to 
minimize the viral load, as well as the drug side-effects and the optimal system is solved numerically to 
determine optimal efficacies of the drugs. Chakrabarty ll20l extended the results in |[T9l by considering a 
clinically validated functional form for the interferon efficacy and hence determined the optimal efficacy 
of ribavirin. Martin et al. |[2TI in a recent paper examine a three compartment model for HCV, involving 
the susceptible, chronically infected and treated injecting drug users (IDUs). They determine an optimal 
treatment programme over a 10 year period taking into account several biomedical and economic objectives. 

2 Model for Monotreatment with IFN 

We first consider a three dimensional model given by Dahari et al. f9\ based on an earlier pioneering 
model of Neumann et al. 17j|. The model describes the dynamics of the the uninfected and infected hepato- 
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cytes as well as the HCV and is expressed as a system of thi^ee coupled ordinaiy-differential equations: 
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The uninfected hepatocytes (T) are being produced at a rate s and proliferate logistically at a rate r, ac- 
companied by a natural death rate of d. The homeostatic liver mechanism is incorporated through Tmax. 
which is the maximum hepatocyte count in the liver. HCV {V) is assumed to infect the hepatocytes at a 
rate /?, thereby producing infected hepatocytes (/). The infected hepatocytes are also assumed to proliferate 
logistically at a rate r and have a natural death rate of 6. In the absence of any kind of treatment, the infected 
hepatocytes produce HCV at a rate p, which has a clearance rate of c. The production of HCV is lowered 
by a factor of (1 — e^) upon the administration of IFN treatment, where is the efficacy of IFN. The model 
excludes the role of IFN in blocking the infection, as it was observed to have minimal effect in this case 
lUlIll- This model admits two steady states 191, namely the uninfected state. 
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Under the physiological conditions r > d and s < dTmax> it can be shown ||9l that there is a transcritical 
bifurcation at, 

— i^m^. — 

In other words, the uninfected steady state is stable if (1 — gp) < Cj*, while the infected steady state is stable 
if(l-ep) >Ci*. 



3 Optimal Control Problem for Monotreatment Model 

In this section, we present an optimal control problem motivated by biomedical considerations. The goal 
is to formulate the problem such that, the integral and terminal values of HCV (V), as well as that of the 
infected hepatocytes (/) are minimized with respect to the control (e^) (on the lines of ifTSl ). The problem 
also incorporates the necessity of minimizing the therapeutic side effects of IFN. Keeping these goals in 
mind, the objective functional, to be minimized, is defined as, 

J = \ [S22l(.tff + SssVitff] + ]- r [Q22l(.tf + QssVitf + Rep{tf] dt, 

J ti 
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where, 5^ and Qu represent the cost coefficient associated with the respective variables. For the purpose of 
notational convenience, we define the state variables as xi = T,X2 = I ,x^ = V and the control variable as 
u = Ep. The state equations in the new variables are given by, 

, /, Xi+X2\ 

xi = s + rxi I 1 — 1 — axi — pxix^ 

\ ^ max / 

, Xi+X2\ 

X2 = PX1X3 + rx2 1 - 0x2 



max 



±3 = {1 - u)pX2 - CX3. (3.1) 

By defining the State vector X = (xi X2 xa)^, the above equations can be written as x = f(x(t), t). 
The initial condition, x(tj) = xq, is taken to be the infected steady state concentrations before the initiation 
of the treatment (when = u = 0) ||9l. 



The above problem can be viewed as a part of a more general setting as described below 
To find an admissible optimal control u* : [ti,tf] [u,u], which for the system x(t) = f (x(i), u(i), i), x(tj) 
xo, minimizes the cost functional 

J = ^x^(t^)5x(t/) + i [xT(t)gx(t) + u^(t)i?u(t)" 

J ti 

= ;C(x(t/),t/) + / C{:>^{t),u{t),t)dt. 

The Hamiltonian, defined by using the dynamic constraint f(x(t), u{t),t) and the Lagrangian i2(x(t), u(t), t) 
through the adjoint vector A-^ (with the same dimension as the state vector), is as follows, 

'H{^{t),u{t), X{t),t) = C{x{t),u{t),t) + A^(t)f (x(t), u(t), t). 

Using Pontryagin's minimum principle, the necessary conditions (in terms of the Hamiltonian) for u* to be 
an optimal control are 



(a) = VA?^(x*(t),u*(t),A*(t),t) =f(x*(^),^x*(^),^)■ 

(b) X*(ti) =xo. 

(c) A*(t) = -[Vx?^(x*(t),u*(0,A*(t),0]'" = -[Vx/:(x*(t),u*(0,i)]^-[Vxf(x*(0,u*(t),t)]"^A(i). 

(d) X*itf) = V,/C i^*itf),tf) = S^*itf). 

(e) Vu'H{^*it),u*{t),X*{t),t) = Vu/:(x*(t),u*(t),t) + A(t)TVuf(x*(i),u*(t),t) = 0. 



The optimal control system thus, comprises of a coupled forward state equation and a backward adjoint 
equation, along with the regular control. This problem, being nonlinear and coupled in nature, needs to be 
solved using concun^ent and iterative numerical procedures. In this paper, the solution is approximated by 
the steepest gradient method |[T8ll22l . 

1. Approximation starts with an initial guess for the control u{t) = u^^^ (t). 

2. Using this control u('^)(t), x^'^) (t) is found by numerically solving x = f (x(t), u^'^) (i), i), x(tj) = xq, 
which is then used to solve A = -Vx'H(x(°) (t), u(°) (t), X{t),t), X{tf) = 5xW(t/) to obtain X^^\t). 
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3. The gradient Vu^(x''°^ (t), u'^''^ (t), A^*^^ it),t) is evaluated to obtain the updated conti"ol as follows, 
uW(t) = uW(t) -r(0)Vu?^(x(°)(t),u(0)(t),A(°)(t),t) 

The above iterative process for the /c + 1-th iteration is given by, 

u('=+i)(t) = uW(t) - r (A;) Vu?^(xW(i),u('=Hi),A('=) 

This procedure is repeated until some convergence criterion, such as | ju^'^"^^) (t) — u^'^^ {t)\\ < tol or J'^''^ — 
j{k+i) ^ jqJ^ satisfied (the latter was used in this paper). Runge-Kutta method of order 4 is used for 
numerically integrating x and A. 



Parameter 


Value (g;] 


s 


1.0 cell ml-i day^i 


d 


0.01 day-i 


P 


2.9 virions day^^ 


P 


2.25 X 10^^ ml day~^ virions"^ 


c 


6.0 day"^ 


5 


1.0 day-i 


r 


2.0 day-i 


T 

max 


3.6 X lO'^ cells ml-i 



Table 1 : Parameter values for the optimal control simulation 

4 Results and Discussion for Monotreatment Model 

We implemented the above procedure using MATLAB ™. The optimal control is approximated using 
the steepest gradient method for the parameter values given in Table [T]f9^]. The efficacy Ep can theoretically 
lie between and 1. Accordingly, u is set to be 0, which corresponds to zero effectiveness of the drug. The 
value of It, however, is taken to be less than 1, for biomedical observations suggest that a perfect efficacy is 
unlikely to be achieved. A more mathematical justification can be presented using pharmacokinetic models 
(such as Powers et al. |[23l ). where the term for efficacy cannot become 1 for all possible values of drug 
doses. 

The optimal control (u*) is simulated over a period of 50 days and is presented in Figure [T] The cost 
coefficients used were S22 = S33 = 10~^Q22 = Q33 = 10"^ and R = 0.1. It can be observed from the 
figure that, the optimal IFN efficacy remains at n = 0.98 for about three weeks and gradually decreases to 
about 0.7 over the next 4 weeks. The same simulation, run over a period of ten weeks, exhibits a slower 
rate of decline in the optimal efficacy. However, the time window for which the optimal efficacy remains 
at u, is roughly unchanged. The simulations are repeated for various values of matrices S, Q and R, for 
various treatment time windows. The results indicate that in responsive patients, the length of this period 
of administering high dosage does not vary by much as long as the cost coefficient R is relatively small. 
Higher values of R result (as one expects) in a shorter period of high dosage. This is of crucial importance 
when deciding on an optimal therapy of drugs with minor side-effects. 

The dynamics of the uninfected (T) and infected (I) hepatocyte concentrations and HCV (V) levels 
under the application of optimal efficacy (Figure [D are given in Figure |2] The viral load shows a triphasic 
decline, which is in line with the biomedical observations ||9l. The optimal treatment exhibits a rapid first 
phase decline in the viral load during the first couple of days, followed by a stable concentration over the 
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next couple of weeks and finally a faster decline for the rest of the treatment period. These observations 
are consistent with the results in |[71|8l|9l. The infected hepatocyte concentration decreases rapidly from the 
third week of the treatment, whereas, the uninfected hepatocyte concentration increases over the course of 
the first three weeks of therapy and remains close to Tmax (maximum liver volume) for the rest of the period. 
Interestingly, the level of optimal efficacy remains at u until a point (three weeks in this case) where the 
viral load and the infected hepatocyte count starts showing a rapid decline, accompanied by the uninfected 
hepatocyte level approaching the maximum liver volume. This is in line with the notion of administering 
high dosage early in the treatment ifTSl and that the dosage should decrease once the results are visible. 
Usually the detection level of virions for HCV patients is taken to be of the order ~ 10^ HCV RNA copies 
per ml IH. We define to be the time at which the virion levels become undetectable. As seen in Figure |2l 
the viral load drops below this level for u* after about 6 weeks (i.e, = 6 weeks). We then examine the state 
dynamics of a large number of sets of sample parameter values under this optimal control u*. The sample 
parameter values are generated using the Latin Hypercube Sampling technique. We implement this sampling 
method to generate a large number of sample points for all parameter values from a multivariate distribution 
(normal, in our case, with a specified mean vector and covariance matrix S) 11241 . For the theoretical and 
computational aspects of this method, the interested reader may refer to the book by Glasserman |[25l . We 

take the mean vector fJ- to be the values specified in Table[Il i.e., fi = {s, d, p, /3, c, 6, r, Tmax)- The diagonal 
covaiiance matrix S is taken as diag (s^, d^, p^, , , 5^ , , T^^^,) ll24ll . Once the sets of parameter values 
are generated they are tested for positivity and physiological conditions (r > d and s < dT^s.x 10) till 
exactly 500 such sets are accepted. The sample sets, which failed this test, are rejected. We, then, examine 
the dynamics of the state variables (by numerically approximating Equation (13.1b ) for all the accepted sample 
parameter sets under the application of the optimal therapy u*. We keep track of the time (t^) required for 
the viral load to fall below detection level of 10^ for each sample set. The cumulative distribution of td 
is shown in Figure [3l The figure shows that, about half of the cases exhibit a short-term response to the 
optimal therapy, which obviously includes the acute cases, where the virus clears out in^espective of whether 
any treatment is administered or not ll24l . 

We also perform Spearman's test to check for monotonic relations between each of the sample parameters 
sets and the time required {tj,) for the viral load to drop below detection level. This test indicates a strong 
negative correlation between td and the death rate of infected hepatocytes 6. A reason for this could be the 
strong negative correlation between 6 and the transcritical bifurcation (C^) over the same distribution |[24l . 
Biologically, the cases, which present with a higher death rate of infected hepatocytes, ai^e more likely (as 
compared to other parameters) to respond to therapy. Also, td shows a mild monotonic correlation with r. A 
higher rate of proliferation results in higher infected hepatocyte levels, which in turn supports the production 
of virions. Therefore, the cases which present with a high rate of proliferation might not achieve short-term 
response. 



In this section, we present an extended model which includes the therapeutic effects of ribavirin to the 
model discussed in Section 2. This model, for combination treatment with IFN and ribavirin, comprises of 
four coupled ODEs: 



5 Model for Combination Treatment with IFN and Ribavirin 
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Ribavirin works by rendering a fraction of the newly produced virions non-infectious. The virion population 
is divided into two different virion populations, namely infectious {Vi) and non-infectious {Vni), as a result 
of administration of ribavirin with an efficacy of p. This model also admits two steady states namely the 
uninfected state, 



and the infected state, 

rp{i) 



2r 



r-d+ \ {r-df + 



4rs 



NI 



0. 



4:STrr 



.^2 



B 



(,) _ (l-ep)(l-p)p/» pV} 



where. 



V, 



^ ^ (l-ep)(l-p)p/3Tniax ^ 
cr ' 



dB 



+ 1 



Under the physiological conditions r > d and s < dT^a.^, it can be shown 191 that there is a transcritical 
bifurcation at, 

- - = — ^pfi^ — = 

In other words, the uninfected steady state is stable if (1 — ep)(l — p) < C| while the infected steady state 
is stable if (1 - ep)(l - p) > C|. 



5. 1 Optimal Control Problem for Combination Treatment with IFN and Ribavirin 

We define a slightly modified objective functional for the four equation model outlined in the previous 
section. In defining the objective functional, we consider the integral values of the infected hepatocyte 
concentration (/) and infectious viral load (V/), as well as the terminal values of infected hepatocytes and 
both the virion populations (Vj and Vni)- The integral value of the non-infectious viral load is not taken into 
account {i.e., has zero cost coefficient), since it does not contribute towards the infectivity of hepatocytes 
during the course of the treatment. Also, the functional incorporates the need to minimize the therapeutic 
side-effects of both IFN and ribavirin. Under these biomedical considerations, the objective functional for 
combination therapy is defined as, 

J = l[S22litff + Ss3Vl{tff + SuVNl{tff] 

+ \ P [Q22l{tf + QssViitf + Rnepitf + R22p{tf] dt, 

The state vector for this model is defined as x = (xi X2 X3 ^4)^, where xi = T,X2 = I, a^s = Vj and 
X4 = Vnj. Similarly, the control vector is defined as u = (ui ^2)^, where ui = ep and U2 = p- The 
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system described in Equation (15.11) is written as x = /(x(t), u{t),t) with the initial condition x(tj) = xq, 
which, as in case of the first model, is taken to be the pre-treatment (ui = U2 = 0) infected steady state 
concentrations. The subsequent mathematical and numerical treatment of the problem is analogous to the 
approach adopted for the monotreatment model. 

5.2 Results and Discussions for Combination Treatment with IFN and Ribavirin 

As in the case of monotreatment, the parameter values are taken from Table [T] The optimal control 
problem is solved using the steepest gradient method described in Section 2. The cost coefficients used for 
numerical implementation ai^e S22 = S-^s = = 10~^, Q22 = Q33 = 10~^ and Ru = R22 = 0.1. 

The optimal combination therapy u* and the dynamics of the state variables under u* are presented 
in Figure (HJl and ^ respectively. As it can be seen from the figures, the patient response, as a result of 
combination therapy, is better as compared to monotreatment regimen. 

It is observed from Figure ^ that, the viral load drops below the detection level (10^) after about 30 days 
of combination treatment as compared to 40 days in case of monotreatment. In addition, decline in infected 
hepatocyte levels in this case is steeper, in comparison with the decline under monotherapy. The optimal 
control remains at u = (0.98, 0.98)^ for about three weeks and then gradually decreases over the rest of 
the period, as it was seen in case of monotherapy as well. If the treatment window is varied, length of this 
period of high efficacy remains roughly unchanged for smaller values of R, which reiterates the point made 
earlier that, in case of drugs with minor side effects, the duration of high dosage remains approximately the 
same. As was the case with monotherapy, the infectious viral load (Vj) shows a triphasic decline during 
the period of administering high dosage, while the uninfected hepatocytes (T) exhibit an increase in levels. 
Beyond this period, the uninfected hepatocyte count remains close to Tmax- The infected hepatocyte count 
(/) shows a concurrent constant level followed by a decline. This observation is consistent with the idea 
that, high doses are administered until visible physical progress is achieved, after which, dosage is reduced. 

We use the Latin hypercube sampling technique described in Section 3 to generate sample parameters, 
using the same ^ and S as before along with positivity and physiological restrictions ||9l for this model 
as conditions for acceptability. The cumulative distribution for is determined for the accepted sample 
parameter sets and is shown in Figure (O. The graph for this distribution shows an improvement in patient 
response (60% after 50 days) under combination therapy, as compared to monotreatment success rate of 
about 50%. Spearman's tests, between sample parameters and td, show that a strong negative correlation ex- 
ists between S and t^. This is consistent with our earlier work |[24l . where 5 was found to be the most crucial 
parameter in determining the patient response. Also, a moderate positive coiTclation was observed between 
the rate of proliferation r and t^. A higher proliferation of infected hepatocytes supports the replication of 
HCV, a result also noted in |[24l . 

6 Conclusion 

In this paper, we consider two well established models for HCV dynamics, incorporating the efficacy 
of the drugs administered. The effects of the drugs as manifested by this model, vary greatly amongst 
individual patents, depending upon the parameter variation. A transcritical bifurcation condition which 
identifies the stability criterion for the uninfected and infected steady state is presented. Optimal treatment 
protocol is determined for both the models, one for monotreatment with IFN and the other for a combination 
therapy of IFN and ribavirin. The optimal therapy shows an initial high level followed by a decline once 
the hepatocyte gradually gets restored. The results, obtained in this paper, are consistent with biomedical 
observations relating therapeutic protocols, especially, the notion of administering high dosage in the early 



8 



stages of infection and gradual lowering of the dosage after the patient has shown signs of recovery. 

With the goal of examining the consequences of administering the optimal efficacy of the two drugs 
in case of several patient scenarios, we generate a large set of patient parameters using Latin hypercube 
sampling. The state dynamics for the generated parameter sets are studied under the optimal therapy. The 
results indicated a better response in case of combination therapy as compared to monotherapy. Statistical 
tests, presented in this paper, show that, patient recovery is highly influenced by the death rate of infected 
hepatocytes. Also a higher proliferation rate abets the viral replication. 

Such statistical analysis can be helpful in determining the optimal therapy on an individual basis, by 
recognizing the major factors which influence patient response. We believe that the analysis presented in 
this paper, combined with pharmacokinetics studies, could play an important role in developing improved 
HCV treatment regimen. 
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